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1 Introduction 



The numerical simulation of quantum field theories with fermions is an interesting and 
difficult computational problem. The basic difficulty is the necessary calculation of very 
large determinants, the determinants of fermion matrices, which can only be achieved by 
some stochastic procedure with the help of auxiliary bosonic "pseudofermion" fields. 

In hybrid Monte Carlo algorithms the number of pseudofermion fields corresponds to 
the number of fermion field components. The evolution of the pseudofermion fields in 
the updating process is realized by discretized molecular dynamics equations. The error 
implied by the finite length of discretization steps is corrected for by a global accept-reject 
decision which involves a fermion matrix inversion. The ingredients of the two-step multi- 
bosonic algorithm [|IJ, which will be considered in this contribution, are somewhat different 
but still in a general sense similar. The number of pseudofermion fields is multiplied by 
the order of a polynomial approximation of some negative power x~ a of the fermion 
matrix. These auxiliary bosonic fields are updated acording to a multi-bosonic action || 
by using simple methods known from bosonic quantum field theories, as heatbath and 
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overrelaxation. The error of the polynomial approximation is corrected also here in a 
global accept-reject decision. This is realized in a so called noisy correction step || by 
using a better polynomial approximation, which realizes a kind of generalized inversion 
of the fermion matrix. 

In my talk I review some recent developments of the multi-bosonic algorithms. The 
performance of the two-step multi-bosonic algorithm is illustrated in a recent large scale 
Monte Carlo simulation of the supersymmetric Yang- Mills theory 0, ||, which is being 
performed at NIC, Jiilich. This application shows that the algorithm is able to cope with 
the difficulties arising at nearly zero gluino masse in reasonably large physical volumes. 



2 Multi-bosonic actions 

The multi-bosonic representation of the fermion determinant @ is based on the approxi- 
mation 

where Nf (Dirac-) fermion flavours are considered and the polynomial P n satisfies 

lim PJx) = x~ N t /2 (2) 

in an interval x G [e, A] . This interval is chosen such that it covers the spectrum of the 
squared hermitian fermion matrix Q 2 . The hermitian matrix Q is defined from the original 
fermion matrix Q 

Q = l5 Q = & (3) 
which is assumed here to satisfy the relation 

Q f = 75Q75 ■ (4) 

Note that in ([]]) only the absolute value of the determinant is taken which leaves out its 
sign (or phase). This can be taken into account at the evaluation of expectation values 
by reweighting. 

For the multi-bosonic representation of the determinant one uses the roots of the 
polynomial rj, (j = 1, . . . , n) 

n 

Pn(Q f Q) = Pn(Q 2 ) = r l[(Q 2 - rj) . (5) 

i=i 

Assuming that the roots occur in complex conjugate pairs, one can introduce the equiv- 
alent forms 

n n 

Pn(Q 2 ) = roUKQ ± h) 2 + u 2 ) = r U(Q - p*)(Q - Ps ) (6) 
0=1 3=1 



where Tj = (pj + ivj) 2 and pj = pj + ivj. With the help of complex scalar (pseudofermion) 
fields Qj X one can write 

n 

det P n (Q t Q)- 1 a J] det[(Q - p*)(Q - p,)]" 1 

oc / exp { - f: E ®ty HQ ~ Pj)(Q ~ Pi)]* ^ 1 • ( 7 ) 
{ 3 I ^ J 

The exponent in fl7|) is the (negative) multi-bosonic action. Since it is quadratic in Q, its 
locality properties are inherited from the fermion matrix Q. For instance, if Q has only 
nearest neighbour interactions then the multi-bosonic action (ffi) extends up to next-to- 
nearest neighbours. 

The multi-bosonic representation of the fermion determinant (0) can be used for a 
Monte Carlo procedure in terms of the pseudofermion fields &j X . The difficulty for small 
fermion masses in large physical volumes is that the condition number A/e becomes very 
large (10 4 — 10 6 ) and very high orders n = O(10 3 ) are needed for a good approximation. 
This requires large storage and the autocorrelation becomes bad since it is proportional 
to n. An additional question is how to control the systematic errors introduced by the 
polynomial approximation in ([I]). In principle one has to perform the limit to infinite 
order of the approximation n — > oo, which means in practical terms a need to investigate 
the dependence of the results on n. 

Several solutions for eliminating the systematical errors due to the finite order of 
approximation n in ([!]) are possible. For instance, one can calculate a correction factor 
from the eigenvalues of Q 2 and introduce a corresponding global Matropolis accept-reject 
step in the updating |J. The necessary calculation of the eigenvalues leads, however, to 
an algorithm growing with the square of the number of lattice points. 

A better solution is to apply a noisy correction step |§ which is especially simple in 
case of Nf = 2 flavours when an iterative inversion is sufficient M. This can be generalized 
to an arbitrary number of flavours in the two-step multi-bosonic scheme [p]] which will be 
discussed in detail in the next section. The special case of Nf = 1 flavours can be dealt 
with in a non-hermitian version [[7], |S| applied directly to Q, insted of Q^Q in (TJ). This 
works well for heavy fermion masses when the spectrum of Q can be covered by en 
ellipse but it would be very cumbersome for small fermion masses where the eigenvalues 
are surrounding zero. 

Another possibility to perform the corrections of the systematic errors of the poly- 
nomial approximation is reweighting in the expectation values. For the special case of 
Nf = 2 this has been introduced in the polynomial hybrid Monte Carlo scheme [ lOf . The 
case of arbitrary Nf can be solved by an appropriate polynomial approximation, which 
can also be implemented in the two-step multi-bosonic approach (see section [3.2|) . 



Other approaches for eliminating the systematic errors are possible, for instance, by 
choosing some specific ways of optimizing the approximate polynomials (see |TIl 0). It is 
also possible to combine the multi-bosonic idea with other methods of dynamical fermion 
simulations. To review all attempts would take too much time and most of the proposals 
were not yet tested in really large scale simulations. In what follows I shall concentrate 
on the two-step multi-bosonic scheme which proved to be efficient in recent simulations 
of SU(2) Yang-Mills theory with light gluinos Qg g. 



3 Two-step multi-bosonic algorithm 

The dynamical fermion algorithm directly using the multi-bosonic representation in fl7|) 
has difficulties with large storage requirements and long autocorrelations. One can achieve 
substantial improvements on both these problems by introducing a two-step polynomial 
approximation Jl|, H|. In this two-step approximation scheme ([|) is replaced by 

^pe>{x)P£\x) = x^' 2 , x G [e, A] . (8) 

The multi-bosonic representation is only used for the first polynomial PfM wich provides 
a first crude approximation and hence the order n\ can remain relatively low. The correc- 
tion factor P^ is realized in a stochastic noisy correction step with a global accept-reject 



decision during the updating process (see section |3.1|) . In order to obtain an exact algo- 
rithm one has to consider in this case the limit ^2^00. For very small fermion masses 
it turned out more practicable to fix some large n-i and perform another small correction 
in the evaluation of expectation values by reweighting with a still finer polynomial (see 
section 



3.1 Update correction: global accept-reject 

The idea to use a stochastic correction step in the updating ||, instead of taking very 
large polynomial orders n, was proposed in the case of Nf = 2 flavours in H]. Nf = 2 is 
special because the function to be approximated is just x~ l and P$(x) can be replaced 
by the calculation of the inverse of xP^(x). For general Nf one can take the two-step 
approximation scheme introduced in ]T[ where the two-step multi-bosonic algorithm is 
described in detail. The theory of the necessary optimized polynomials is summarized in 
section § following (13| . 



In the two-step approximation scheme for Nf flavours of fermions the absolute value 
of the determinant is represented as 



|det(Q)|^ ~ j—-^- . (9) 

det P«(Q 2 ) det Pi?(Q^ 
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The multi-bosonic updating with ri\ scalar pseudofermion fields is performed by heatbath 
and overrelaxation sweeps for the scalar fields and Metropolis sweeps for the gauge field. 
After a Metropolis sweep for the gauge field a global accept-reject step is introduced in 
order to reach the distribution of gauge field variables [U] corresponding to the right hand 
side of (|9[). The idea of the noisy correction is to generate a random vector 77 according 
to the normalized Gaussian distribution 



and to accept the change [U'] <— [U] with probability 



(10) 



mm{l,A( V ;[U'}^ [£/])} , (11) 

where 

A{ m [u'\ - [U]) = ex P {-^p£\Q[uT)v - rfP£(Qiu] 2 )v} ■ (12) 

The Gaussian noise vector rj can be obtained from rj' distributed according to the 
simple Gaussian distribution 

e-i/V 

(13) 



by setting it equal to 

v = Pg(Q[V\ 2 )-fy . (14) 

In order to obtain the inverse square root on the right hand side of (|H|), we can proceed 
with polynomial approximations in two different ways. The first possibility was proposed 
in [|l| with x = Q 2 as 

Pg(x)-i * R n3 (x) ~ x N ^S ns [P£(x)\ . (15) 

Here 

S„ s (P)~P^ (16) 

is an approximation of the function P^ on the interval P G [\~ Nf / 2 , e~ N ^ 2 \. The poly- 
nomial approximations R m and S ns can be determined by the same general procedure as 
pW anc i it turns out that these approximations are "easier" in the sense that for a 
given order higher precisions can be achieved than, say, for P$. 

Another possibility to obtain a suitable approximation for (|14|) is to use the second 
decomposition in (|6]) and define 

n.2 

^: /2) (Q)-v^n(Q-p,) , p${q 2 )=p£J 2 Xq)'p% 2 \q) ■ (17) 
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Using this form, the noise vector 77 necessary in the noisy correction step can be generated 
from the gaussian vector rj according to 

v = p£/ 2) {Q)-W, (is) 

where P^ 2 \Q)~ l can be obtained as 



P^iQY 1 = - p ns(Q 2 ) p nl /2} (QV • (19) 



In the last step P n3 denotes a polynomial approximation for the inverse of P$ on the 
interval [e, A]. Note that this last approximation can also be replaced by an iterative 
inversion of P^(Q 2 ). However, tests showed that the inversion by a least-squares opti- 
mized polynomial approximation is much faster because, for a given precision, less matrix 
multiplications have to be performed. 

In the simulation with light dynamical gluinos 0, || mainly the second form in (|T^ ) - ([T9|) 
has been used. The first form could, however, be used as well. In fact, for very high orders 
ri2 or on a 32-bit computer the first scheme is better from the point of view of rounding 
errors. The reason is that in the second scheme for the evaluation of P^J {Q) we have to 
use the product form in terms of the roots pj in fllTD- Even using the optimized ordering 
of roots defined in Jl[ |l3j, this is numerically less stable than the recursive evaluation 
according to (|25l), (|3ll). If one uses the first scheme both P® 1 in ( |12|) and R n3 in ([14]) - 



fllq) can be evaluated recursively. Nevertheless, on a 64-bit machine both methods work 
well and in case of (|19"D the determination of the least-squares optimized polynomials is 
somewhat simpler. 

The global accept-reject step for the gauge field can be performed after full sweeps 
over the gauge field links. A good choice for the order m of the first polynomial PW is 
such that the average acceptance probability of the noisy correction be near 90%. One 
can decrease n\ and/or increase the acceptance probability by updating only some subsets 
of the links before the accept-reject step. This might be useful on lattices larger than the 
largest lattice 12 3 • 24 considered in ||. 



3.2 Measurement correction: reweighting 

The multi-bosonic algorithms become exact only in the limit of infinitely high polynomial 
orders: n — > 00 in (Q) or, in the two-step approximation scheme, n<i — > 00 in Instead of 
investigating the dependence on the polynomial order by performing several simulations, 
it is practically better to fix some high order for the simulation and perform another 
correction in the "measurement" of expectation values by still finer polynomials. This is 
done by reweighting the configurations in the measurement of different quantities. In case 

6 



of Nf = 2 flavours this kind of reweighting has been used in JT0| within the polynomial 
hybrid Monte Carlo scheme. As remarked above, Nf = 2 is special because the reweighting 
can be performed by an iterative inversion. The general case can, however, also be treated 
by a further polynomial approximation. 

The measurement correction for general Nf has been introduced in WM. It is based 
on a polynomial approximation which satisfies 

n limPW(x)^(x)P|J(x) = x~ N f' 2 , x G [e', A] . (20) 

The interval [e f , A] can be chosen, for instance, such that e' = 0, A = A max , where X m ax is 
an absolute upper bound of the eigenvalues of Q^Q = Q 2 . In this case the limit n 4 — > oo 
is exact on an arbitrary gauge configuration. For the evaluation of Pjg one can use n^- 
independent recursive relations (see section f|), which can be stopped by observing the 
convergence of the result. After reweighting the expectation value of a quantity A is given 
by 



where rj is a simple Gaussian noise like rf in ([□]). Here (. . .)u iV denotes an expectation 
value on the gauge field sequence, which is obtained in the two-step process described in 
the previous subsection, and on a sequence of independent 77's. The expectation value 
with respect to the 77-sequence can be considered as a Monte Carlo updating process with 
the trivial action S v = rfrj. The length of the ^-sequence on a fixed gauge configuration 
can be, in principle, arbitrarily chosen. In praxis it can be optimized for obtaining the 
smallest possible errors. 

The application of the measurement correction is most important for quantities which 
are sensitive for small eigenvalues of the fermion matrix Q^Q. The polynomial approxi- 
mations are worst near x = where the function x~ N f' 2 diverges. In the exact effective 
gauge action, including the fermion determinant, the configuration with a small eigen- 
value A are suppressed by A. N ^ 2 . The polynomials at finite order are not able to provide 
such a strong suppression, therefore in the updating sequence of the gauge fields there 
are more configurations with small eigenvalues than needed. The exceptional configura- 
tions with exceptionally small eigenvalues have to be supressed by the reweighting. This 
can be achieved by choosing e' = and a high enough order n^. It is also possible to 
take some non-zero e' and determine the eigenvalues below it exactly. Each eigenvalue 
A < e' is taken into account by an additional reweighting factor A N f^ 2 P^(A)Pj 2 \A). 
The stochastic correction in (pij ) is then restricted to the subspace orthogonal to these 
eigenvectors. Instead of e' > one can also keep e' = and project out a fixed number of 
smallest eigenvalues. 
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Let us note that, in principle, it would be enough to perform just a single kind of 
correction. But to omit the reweighting does not pay because it is much more comfortable 
to investigate the (small) effects of different n 4 values on the expectation values than to 
perform several simulations with increasing values of 7i2- Without the updating correction 
the whole correction could be done by reweighting in the measurements. However, in 
practice this would not work either. The reason is that a first polynomial with relatively 
low order does not sufficiently suppress the exceptional configurations. As a consequence, 
the reweighting factors would become too small and would reduce the effective statistics 
considerably. In addition, the very small eigenvalues are changing slowly in the update 
and this would imply longer autocorrelations. 

A moderate surplus of gauge configurations with small eigenvalues may, however, 
be advantageous because it allows for an easier tunneling among sectors with different 
topological charges. For small fermion masses on large physical volumes this is expected 
to be more important than the prize one has to pay for it by reweighting, provided that 
the reweighting has only a moderate effect. 



4 Least-squares optimized polynomials 

The basic ingredient of multi-bosonic fermion algorithms is the construction of the nec- 
essary optimized polynomial approximations. The least-squares optimization provides a 
general and flexible framework which is well suited for the requirements of multi-bosonic 



algorithms [13|]. In the first part of this section the necessary basic formulae are col- 
lected. In the second part a simple example is considered: in case of an appropriately 
chosen weight function the least-squares optimized polynomials for the approximation of 
the function x~ a are expressed in terms of Jacobi polynomials. 



4.1 Definition and basic relations 

The general theory of least-squares optimized polynomial approximations can be inferred 
from the literature |16|, [17]]. Here we introduce the basic formulae in the way it has been 
done in fl3| for the specific needs of multi-bosonic fermion algorithms. We shall keep the 



notations there, apart from a few changes which allow for more generality. 

We want to approximate the real function f(x) in the interval x G [e, A] by a polynomial 
P n {x) of degree n. The aim is to minimize the deviation norm 




5 n = {NT* / dxw(x) 2 [f(x)-P n (x)Y\ . (22) 
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Here w(x) is an arbitrary real weight function and the overall normalization factor N €t \ 
can be chosen by convenience, for instance, as 

N e ,\ = dxw(x) 2 f(x) 2 . (23) 

A typical example of functions to be approximated is f(x) = x~ a /P(x) with a > and 
some polynomial P(x). The interval is usually such that < e < A. For optimizing the 
relative deviation one takes a weight function w(x) = 

It turns out useful to introduce orthogonal polynomials (/i = 0, 1, 2, . . .) satis- 

fying 

rX 

J dxw(x) 2 <$> a (x)<$> u (x) = S^Qu . (24) 
and expand the polynomial P n (x) in terms of them: 

n 

P n (x) =Y, d nv®v{x) . (25) 
u=0 

Besides the normalization factor q u let us also introduce, for later purposes, the integrals 
p u and s„ by 

/A rX rX 

dxw{x) 2 Q v {x) 2 , Pu = dxw{x) 2 Q u {x) 2 x , s v = / dxw(x) 2 x u . (26) 

It can be easily shown that the expansion coefficients d nu minimizing 5 n are indepen- 
dent from n and are given by 

dnu = d u = — , (27) 

where 

rX 

bu = j dxw(x) 2 f(x)<5>u(x) . (28) 

The minimal value of S 2 is 

n 

# = 1 " ^a E <*A • (29) 

i;=0 

The above orthogonal polynomials satisfy three-term recurrence relations which are 
very useful for numerical evaluation. The first two of them with fj, — 0, 1 are given by 

$ (aO = l, $i(x)=a;--. (30) 

The higher order polynomials $n(x) for /x = 2, 3, . . . can be obtained from the recurrence 
relation 

<Vi(x) = (x + ^)^(x)+7 M _ 1 $ M _ 1 (a;) , (// = 1, 2, . . .) , (31) 
where the recurrence coefficients are given by 

/?„ = --> 7,-i = — ^ • (32) 
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Using these relations on can set up a recursive scheme for the computation of the 
orthogonal polynomials in terms of the basic integrals s v defined in (|26|). Defining the 
polynomial coefficients f^ u (0 < v < p) by 

*M = fl Uu^~ V (33) 
the above recurrence relations imply the normalization convention 

/ M o=l, (^ = 0,1,2,...), (34) 

and one can easily show that and satisfy 

n fl 
Qfi = J]] ffj 1 u s 2fi-u i Vfi = f[iv { s 2fi+l-u + ffil s 2fi-u ) • (35) 

The coefficients themselves can be calculated from / n = — si/s Q and ([JT]) which gives 

f 11+1,1 = ffj.,1 + Pfi , 

ffi+1,2 = f/j.,2 + ftfif /J.,1 + ln-1 , 
JfM+1,3 = ffj.,3 + ft/if [i,2 + lfi-lffi-1,1 , 

f/i+1,/1 = ffJ.,Li "I" f^iif ji,ii— 1 "I" 1 li— if '/J-— l,/i— 2 ? 

f ii+i,ii+i (3fif[i,fi ~\~ 7 fi— lj ? fi— i,ii— i ■ 

(36) 



The polynomial and recurrence coefficients are recursively determined by (|34])-(|3"6|). The 
expansion coefficients for the optimized polynomial P n (x) can be obtained from (^) and 

b, = J2U» f X dxw(x) 2 f(x)x^ . (37) 
4.2 A simple example: Jacobi polynomials 

The approximation interval [e, A] can be transformed to some standard interval, say, [—1,1] 
by the linear mapping 

2x — A — e £ , , . 1 , , . . , 

t= x _ e , x = |(A-e) + -(A + e) . (38) 

A weight factor (1 + £) p (l — with p, o > — 1 corresponds in the original interval to the 
weight factor 

W M ^2 = ^ _ e y^ x _ x y _ ^ 

Taking, for instance, p = 2a, o = this weight is similar to the one for relative deviation 
from the function f(x) = x~ a , which would be just x 2a . In fact, for e = these are exactly 
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the same and for small e the difference is negligible. The advantage of considering the 
weight factor in ( |39"D is that the corresponding orthogonal polynomials are simply related 
to the Jacobi polynomials f!8 , [19|| , namely 

*™ (x) = (A - e?v\ r ^ + g + y + 1 ) p(M ( 2 "~ A -^ . (40) 

Our normalization convention (|34]) implies that 

„M = (\- AP+^+h/ F(p + - + + - + 1)r(p + <y + v + l ) (a, \ 

The coefficients of the orthogonal polynomials are now given by 
/CM = t(-)-( £ -A)" ( U " ) ^t" +1>r( Cr + 2 ",Tt!! ■ («) 



1 /jy 

oj=0 



V — U I \ u 



T(p + fi - tu + l)T(p + a + 2p+l) 



In particular, we have 

fft* = 1 , A ( r } = — e — (A — . (43) 



The coefficients /3, 7 in the recurrence relation (|31|) can be derived from the known recur- 
rence relations of the Jacobi polynomials: 

^ , ^ , (^ 2 -p 2 )(A-e) 



2 V ' 2(p + a + 2p)(p + a + 2p + 2) ' 

(p,«t) = _(\-e) 2 p(p + p)(a + p)(p + (T + p) . . 

7m_1 1 J (p + cr + 2p- l)(p + a + 2p) 2 (p + a + 2p + l) ' 1 ' 

In order to obtain the expansion coefficients of the least-squares optimized polynomials 
one has to perform the integrals in (fT?|) . As an example, let us consider the function 
f(x) = x~ a when the necessery integrals can be expressed by hypergeometric functions: 

£ dx (x - e) p (X - xYx»- v - a = 

= (A - e) p + ^ A ^ r ^ + |)^+l) F ( a - p + „, (T + 1;p + a + 2;l-^ . (45 ) 
Let us now consider, for simplicity, only the case e = 0, when we obtain 

1,0,*) = ( ]r} i+ P +*+»-a T (P + (T + V+ 1 ) T ( a + t J ') T (p- a + 1 ) T ( (T + L L + 1 ) 

M ~ l ' r(p + a + 2p + i)r(a)r( P + (7-o i + p + 2) ■ 1 ] 



Combined with (P7|) and fl4TD this leads to 

,yO,<r) _ f 1 ^ \ T(p + a + 2p + 2)r(q + p)r(p - a + 1) 
" 1 ' p!r(p + p + l)r(q)r(p + a-a + p + 2) • 1 J 
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These formulae can be used, for instance, for fractional inversion. For the parameters 
p, a the natural choice in this case is p = 2a, a = which corresponds to the optimization 
of the relative deviation from the function f(x) = x~ a . As we have seen in section 
[4.1| , the optimized polynomials are the truncated expansions of x~ a in terms of the Jacobi 
polynomials p( 2a '°) . The Gegenbauer polynomials proposed in [20] for fractional inversion 
correspond to a different choice, namely p = a = a — |. This is because of the relation 

= If^^'^W • (48) 

Note that for the simple case a = 1 we have here the Chebyshev polynomials of second 
kind: C l n (x) = U n (x). 

The special case a = \ is interesting for the numerical evaluation of the zero mass 



lattice action proposed by Neuberger [pifl . In this case, in order to obtain the least- 
squares optimized relative deviation with weight function w{x) = x, the function x~^ 
has to be expanded in the Jacobi polynomials pQ>°\ Note that this is different both 
from the Chebyshev and the Legendre expansions applied in p2]| . The former would 
correspond to take p(~i~\\ the latter to p(°>°). The corresponding weight functions 
would be [x(\ — x)}~^ and 1, respectively. As a consequence of the divergence of the 
weight factor at x = 0, the Chebyshev expansion is not appropriate for an approximation 
in an interval with e = 0. This can be immediately seen from the divergence of dp 2 ' 2 
at a = \ in (|47|). 

The advantage of the Jacobi polynomials appearing in these examples is that they 
are analytically known. The more general least-squares optimized polynomials defined 
in the previous subsection can also be numerically expanded in terms of them. This is 
sometimes more comfortable than the entirely numerical approach. 



5 Outlook 

The two-step multi-bosonic algorithm has been shown to work properly in a recent large 
scale numerical simulation of light dynamical gluinos Since gluinos are Majorana 

fermions the effective number of (Dirac-) fermion flavours is in this case Nf = |. This 
simulation is demanding because it deals with nearly zero and negative fermion masses 
in reasonably large physical volumes. It turned out possible to investigate the expected 
first order phase transition at zero gluino mass. The algorithm was able to cope with the 
metastability of phases near the phase transition. Up to now the investigated physical 
volumes were not very large: in case of the largest (12 3 -24) lattice the product of the spatial 
extension times the square root of the string tension was L^fo ~ 2.4. The experience 
shows, however, that simulations with light gluinos and L^fo ~ 5 or larger would be 
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feasible with reasonable effort. Another important conclusion in || is that the inclusion 
of the sign of the fermion determinant, actually Pfaffian for Majorana fermions, can be 
achieved by determining the spectral flow of the hermitian fermion matrix Q as a function 
of increasing hopping parameter. 

An interesting application of the two-step multi-bosonic algorithm is the numerical 
simulation of QCD. Up to now no serious effort was taken in this direction but first tests 
will be performed soon. A particularly relevant aspect for QCD is the possibility to deal 
with an odd number of quark flavours (Nf). The popular hybrid Monte Carlo algorithm 
can only be applied for even Nf. For non-even Nf one can use finite step-size molucular 
dynamics algorithms like the one in [23], but then the extrapolation to zero step size is 
an additional difficulty 
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